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o 

We present a Langevin approach to describe the steady-state dynamics of one-dimensional 
granular media fluidized by a vibrating bottom plate. We adopt a linear Langevin equation 
' to describe the motion of the center of mass. Within this framework, we derive analytical 

expressions for several macroscopic quantities. We also predict the power spectrum for the 
height of the center of mass. We find good agreement between our theoretical predictions 
and extensive event-driven molecular dynamics simulations. 
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1. Introduction 

Granular materials fluidized by external vibrations have attracted a great deal of interest 
in the field of nonequilibrium statistical physics. These systems show a variety of fascinating 
phenomena (see, for example, ref. 1), such as convection, surface waves, and size segrega- 
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tion. The transition from a fluidized state to a condensed state, in which the particles move 
collectively with the same period as the external vibrations, has been studied both experi- 
mentally 2-4 ) and numerically. 2, 3,s ) 

Scaling laws of the macroscopic properties (e.g. the center of mass (COM)) in a fluidized 
state are important subjects from the viewpoint of a nonequilibrium steady state (NESS), 
and have also been extensively studied using experiments, 3 ' 6 ' 7 - 1 simulations, 3, 8 ~ 10 ) kinetic the- 
or y 6, n-13) an( j bydrodynamic descriptions. 14 ' 15 ) However, no agreement has been reached 
among the various studies for the scaling relationship of the height of the COM in a NESS. 
This discrepancy is not yet understood and the question remains unanswered. 

To resolve this problem, we focused on a one-dimensional column of inelastic hard rod 
(granular) particles bouncing on a sinusoidally vibrating bottom plate. A one-dimensional 
system allows us to study the pure dynamics of the COM in a fluidized state since we can 
remove the effect of convection within granular beds, which becomes dominant in higher 
dimensions. A one-dimensional system has been investigated in detail previously by several 
authors. Numerical simulations 3 ) have clarified the behavior of the transition from a condensed 
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state to a fluidized state and showed a scaling relationship for the height of the COM in the 
fluidized state. The transition is characterized by a parameter X = N(l — r), where N is 
the number of granular particles and r is the restitution coefficient for collisions between 
particles. If X is greater than the critical value Xq, which is about 3, then the transition 
does not occur regardless acceleration of the vibrating bottom plate; if X < Xc, then the 
transition may occur at a certain magnitude of the acceleration. As an explanation of the 
value of Xc — 3, only connection with the result of an analytical study 16 ) on inelastic particles 
colliding with a wall has been pointed out. The scaling relationship for the height of the COM 
has been shown to be characterized by X and typical velocity of the bottom plate. Theoretical 
considerations 11 ^ based on kinetic theory have shown that analytical solutions of a dissipative 
Boltzmann equation for a one-dimensional system agree well with the results of numerical 
simulations. 

The main purpose of our study was to develop a novel theoretical method to clarify the 
macroscopic behavior of fluidized granular media through the dynamics of the COM without 
using a kinetic equation. Our theoretical method has a wide range of applications to both gas 
states, for which kinetic theory can be applied, and dense liquid states, when the COM of the 
granular media evolves with the same period of the bottom plate. Our basic formulation, which 
consists of a Langevin equation of motion of the COM, was derived by focusing on the force 
acting on the granular media due to the bottom plate. Some macroscopic quantities can be 
derived easily from the solution of this equation. Although several phenomenological parame- 
ters must be estimated by numerical simulation, we confirmed that our compact formulation 
agreed with extensive event-driven molecular dynamics simulations self-consistently. 

We considered a column of iV particles of mass m bouncing on a vibrating bottom plate 
subjected to gravity with acceleration g. The motion of the particles was confined to the 
z axis, and the z direction was opposite to the direction of gravity. Since the diameter of 
the particles plays no role in one-dimensional hard sphere systems, we regarded the particles 
as point particles. The bottom plate oscillated sinusoidally with amplitude Aq and angular 
frequency to; hence, the height of the bottom plate zo(t) at time t was 

z (t) = A sin (ut) . (1) 

Collisions between particles were inelastic with a restitution coefficient r. For simplicity, we 
assumed that collisions between the lowest particle and the bottom plate were elastic. 

Two important timescales occur in this system: the oscillation period of the bottom plate 
t(= 2w/lo), and the macroscopic relaxation time T re \ to the stationary state. If r is comparable 
to T re i , the energy supplied by one stroke of the bottom plate is almost dissipated during the 
period, and the particles will be in a condensed state. Such a condensed state is beyond the 
scope of the present study. In this study, we restricted ourselves to the high-frequency case 
t '/i~rel <C 1, in which the system is in a fluidized state. 
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We performed simulations systematically by changing the number of particles N (TV = 
10, 100, 1000), the restitution coefficient r (r = 0.80 ~ 0.9999), and the maximum acceleration 
of the bottom plate T = Aqlo 2 jg (T = 10 ~ 640). The physical quantities were averaged over 
a long period of time: 1.0 x 10 5 r for N = 10 and 100, and 5.0 x 10 4 r for N = 1000. 

2. A Langevin Approach 

To describe the macroscopic properties of the fluidized state, we started from the equation 
of motion of the COM: 

M^- = -Mg + F b , (2) 

where Z{t) is the height of the COM of the particles at time t, M = Nm is the total mass, 
and the right hand side of eq. (2) represents the sum of all external forces acting on the 
granular particle system: the first term is the gravitational force and the second term F b (t) is 
the external force exerted by the vibrating bottom plate at time t. The long-time average of 
Fb(t) must balance with the gravitational force Mg acting on the column of particles. Hence, 
the central part of this study evaluates the deviation of F b {t) from its long-time average: 
SF(t) = F b (t) - Mg. 

To evaluate 5F(t), we focused on the reaction force of F b (t), i.e., the force exerted on 
the bottom plate by the granular fluid. We drew an analogy between the force acting on 
the bottom plate and the force acting on a fine particle exhibiting Brownian motion while 
immersed in a fluid (see, for example, ref. 17); here the bottom plate on one side of the granular 
fluid corresponds to a Brownian particle in a fluid. 18 ) The force on one side of the Brownian 
particle consists of the pressure, the frictional force, and the random force. We assumed that 
the force on the bottom plate also contained these three characteristic forces: the average force 
Mg on the bottom plate corresponding to the average pressure on the Brownian particle, the 
frictional force accompanied by the relative motion of the bottom plate against the granular 
fluid, and the random force. In addition we have to take into account that time dependence of 
the pressure which is originated from the elasticity of the granular fluid. We assumed here two 
kinds of elastic forces which cause time dependence of the pressure. The first is the reaction 
force when the bottom plate excites a sound wave in the granular fluid. The second is the 
elastic force accompanied by macroscopic motion of the granular fluid; this force is significant 
in our system because the granular fluid of finite length has the slowest oscillating mode with 
a macroscopic time scale. 

Finally, we assumed the following form of 5F(t), consisting of systematic forces and a 
random force: 

6F(t) = -k (Z(t) -Z) + f s (t) - fiV(t) + R{t). (3) 
The first term represents the elastic force accompanied by the slowest mode showing the 
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expansion and contraction for the total length of the column of particles. We assumed this 
was proportional to the deviation of the height of the COM Z(t) from its stationary value 
Z, Z(t) — Z, where Z is defined by Z = YniiT—>oo t> J Z(t)dt. The constant k is the elastic 
constant. The second term in eq. (3) is the elastic force f s (t) accompanied by the excitation 
of a sound wave at the bottom plate. The third term is the frictional force, which we assumed 
was proportional to the COM velocity V(t). The constant p, is the frictional constant. The 
last term R(t) is the random force, which is specified later. 

We estimated the elastic force f s (t) on the basis of hydrodynamic sound-wave theory. 19 ) 
In a normal fluid, sound waves propagate according to a relationship between the pressure 
and the velocity of the fluid. Let us denote a small change in the pressure from its equilibrium 
value by p', a typical velocity of the fluid particles in the wave by v, and the velocity of sound 
by c s . If the condition v <C c s is satisfied, we have a relationship p' = pc s v for a traveling 
plane wave, where p is the constant equilibrium density of the fluid. We assumed that this 
relationship was also satisfied in fluidized granular media under the same condition v <C c s . 
Let us introduce here the global temperature T and the thermal velocity c of the granular 
fluid defined in analogy with those of a normal fluid. The global temperature is related to 
the mean square velocity fluctuation by k^T = m((v — {v)n) 2 )n, where ks is the Boltzmann 
constant and the angle brackets (■■■)at denote the average over all particles. The thermal 
velocity c is defined using the global temperature T as c = ^Jk^T '/m. The velocity of sound 
c s is on the order of the thermal velocity c. The density p is on the order of M/(c 2 /g), where 
c 2 /g is the maximum height reached by a particle launched from the bottom plate with a 
thermal velocity c; this maximum height characterizes the length of the column of particles. 
In the vicinity of the bottom plate, v may be approximated by the velocity of the bottom 
plate vo(t) = Aquj cos(ut). Since f s (t) corresponds to p' at the bottom plate, we have 



where a is a numerical factor on the order of 1 that is used as a curve-fit parameter when we 
compare our theoretical predictions with the results of simulations. 

There are some important consequences derived from the condition v <C c s . First, in 
the vicinity of the bottom plate, the condition v <C c s can be written as vq <^ c. Hence, the 
maximum value of vo(t), Aqlo, must be small compared to c: Aqlo <C c. Secondly, the condition 
v < c s suggests (v) 2 N <C c 2 = (v 2 )n - (v) 2 N — (v 2 )n- Hence, 



where Ek is the stationary value of the kinetic energy per particle defined as Ek = 
W J2iLi Y v i- I n our simulation, we estimated the thermal velocity c according to the relation 



The elastic constant k is related to the angular frequency f2 of the slowest mode of the 



f s {t) = aM-A ujcos(iot), 



(4) 
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macroscopic oscillatory motion by £1 =. sj h j J\d . The period t osc of the oscillation is given by 
t osc = 27r/0. The rate of relaxation to the stationary state // due to friction is defined by 
// = fi/M, which is related to the macroscopic relaxation time r re / by r re ; = Since both 
t osc and r re i characterize the macroscopic change that extends to the full length of the column 
of particles, they must be on the same order as the characteristic time taken for a sound wave 
to travel along the total length of the column of particles. This characteristic time can be 
estimated as c/g because the velocity of sound is on the order of c and the total length of the 
column of particles is on the order of c 2 /g. Thus, we assume that Q and // are on the order 
of g/c: 

$7 = Cl—, u! = u-, (6) 

c c 

where Cl and fi are numerical factors on the order of 1 that are determined by curve-fitting 
the results of our simulations. 

We supposed that the properties of the random force R(t) were the same as those acting 
on a fine particle undergoing Brownian motion in a fluid at temperature T. 17 ) Thus, R(t) is 
stationary Gaussian white noise: 

(R(t))=0, (R(t)R(t'))=m-t). (7) 

We assumed that the fluctuation-dissipation theorem was satisfied: 

I = 2Mn'k B T = 2Mg(imc, (8) 



where we used eq. (6) and the relationship c = y/k^T/m to derive the second equality. 

Collecting the above results, we have the following linear Langevin equation for the COM: 

™ = -*(z-D-Sv + ± + ±. (») 

This equation has the same form as the Langevin equation that describes the forced oscillations 
of a fine particle undergoing Brownian motion in a harmonic potential. The solution to this 
equation has the form 

Z{t)-Z = A (sm(ut + 6)+ f G(t - t')^p-dt' 

+F im (t), (10) 



where 



and 



C = , (ii) 

y/(n 2 - u 2 ) 2 + Gu'u;) 2 
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Fig. 1. A test of the theoretical results given by eq. (14). The accelerations T used in the simulations 
were T — 10, 20, 40, 80, 160, 320, and 640, except for the simulations that experienced an "inelastic 
collapse" (an infinite number of collisions in a finite time). 16 ' 20 ) The dashed line corresponds to 
the theoretical prediction a/ui with a — 1.5. 
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respectively. The function G(t) is given by 

G(t) = sin (wot), (13) 

where ujq = (f2 2 — (///2) 2 ) 1 / 2 . The last term Fi n i(t) consists of those that depend on the 
initial conditions and vanish after a sufficient amount of time. Thus, the term is negligible 
when calculating long-time averages of physical quantities in a stationary state. 

Let us consider a fluidized state with the timescale (r /r re ;) 2 <C 1. This can be rewritten 
by substituting r = 2tt/uj and r re \ = c/g as Cj 2 S> 1, where Cj is defined as Cj = uic/g. In this 
limit, we expanded £ in terms of lj~ 2 : 

C = |(l + Or 2 )). (14) 

Figure 1 gives the simulation results for (asa function of Cj with c = (2i?^'/m) 1 / 2 . As shown 
in previous studies, 2 ' 3 ' n ) the relevant parameter that governs the behavior of the system was 
X = (N— 1)(1 — r) rather than N or r; X is the parameter that measures effective dissipation 
in the column of granular particles. The factor N — 1 in the definition of X represents the 
number of dissipative contacts in the system, where —1 is due to the fact that we assumed 
elastic collisions between the lowest particle and the bottom plate. Data points with the same 
T and X coincided in the figure. For simulations with X < 1, shown by circles and squares, the 
master curve was consistent with the theoretical predictions (14) using a = 1.5. For X > 1, 
the numerical factor a weakly depended on X. We also observed deviations from the master 
curve at large Cj for the simulations with small X. These deviations may be attributable an 
insufficient simulation time to obtain the stationary average. 

Some macroscopic quantities can be calculated using the formula given by eq. (10). Let 
us first consider the power injected by the bottom plate Pb, which has been the subject of 
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Fig. 2. A test of the theoretical result given by eq. (16). The accelerations Y used in the simulations 
were Y = 10,20,40,80,160,320, and 640, except the simulations that experienced an inelastic 
collapse. The dashed line corresponds to the theoretical predictions (a/2) (A u>/c) with a = 1.5. 



recent studies. 6 ' 12 ' 13 ' 21,22 - ) Using the force supplied by the bottom plate F b and its velocity 
vo, the power input P b can be defined by P b = FbVo = lim-r^oo ^ J Q T F b (t)vo(t)dt. Substituting 
F b using eq. (2) in the definition of Pf, and integrating by parts twice, we have 

P b = -Muj 2 lim i f Z(t)v (t)dt. (15) 
Then, substituting eq. (10) into eq. (15), we obtain 

= ^(110(0). (16) 

The results obtained by neglecting terms on the order of Cj~ 2 coincided with the scaling 
predicted by kinetic theories: 6, 12 ' 13 ) P b ~ Mg(A u;) 2 /c. Figure 2 shows that the scaling 
relationship (16) with a factor a = 1.5 agrees well with the simulations in the range when 
Aqlo/c < 1. This result is consistent with the condition Aquj/c <C 1 required for the formula 
given by eq. (4) to be valid. 

In previous studies 6 ' 12 ' 13 ) on the basis of kinetic theory, the scaling relationship for the 
granular temperature T ~ (Aquj) 2 / X was derived from a balance between the power input 
from the bottom plate P b and the rate of energy dissipation due to inelastic collisions between 
particles. If we assume the scaling relationship is satisfied in the parameter range where our 
simulations were performed, we have c ~ A§ujj\/X~. Then, the condition Aquj/c < 1 can be 
rewritten as Vx < 1. This result is consistent with the observation that the simulation data 
agreed well with the theoretical predictions in Fig. 1 if X < 1. 

Next, let us consider the power spectrum Icm for the height of the COM. According to 
the Wiener-Khinchin theorem, this can be calculated analytically from the Fourier transform 
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of the two-time correlation function V'CM(i) denned by 

Y>CM(i) = lim i [ T (5Z(t')5Z(t> + t)) dt' , (17) 

T^oo 1 J Q 

where SZ(t) = Z(t) — Z and the brackets (• • • ) indicate an average over the random force 
R{t). Substituting eq. (10) into eq. (17) and performing the Fourier transform gives 

+ r# — , (18) 

(ft 2 -co' ) 2 + (fico') 2 

where to' is the angular frequency to' scaled by g/c: to' = to'c/g. The first term in eq. (18) gives 
the delta functional peaks at Co' = ±u) while the second term gives a continuous spectrum; 
these characteristics of the power spectrum have already been observed in refs. 2 ' 3 ' 6 ) In Fig. 3, 
the scaled power spectrum successfully collapsed onto a single master curve that agreed well 
with the curve given by the second term in eq. (18) having the curve-fit numerical factors 
ft = 2.0 and Cl = 1.5. We stress here that the curve given by the second term in eq. (18) is 
independent of the system parameters N, r, and T. 

3. Conclusions 

This paper has studied a fluidized state of a one-dimensional vibrated granular media, 
using a Langevin approach. We derived a Langevin equation of motion of the COM assum- 
ing the forces acting on the granular fluid due to the bottom plate under phenomenological 
consideration. Using the solution of the Langevin equation, we were able to obtain analytical 
expressions for several quantities, the amplitude of the motion of the COM, the power input 
from the bottom plate, and the power spectrum for the height of the COM. 

In our theory, properties of the granular fluid are characterized only by two macroscopic 
quantities, the total mass M and the thermal velocity c. Because inelastic collisions between 
particles are not explicitly handled in our theory, the restitution coefficient r and the number 
of dissipative contacts N — 1 in the system do not appear in the Langevin equation (9). 
These parameters control the thermal velocity c and influence indirectly other macroscopic 
quantities. In the previous theoretical studies, 6,12 ' 13 ) the scaling relationship of the granular 
temperature T, which is related to the thermal velocity c by c = -y/fceT/m, on the system 
parameters including r and N has been obtained in the following way: They estimated on 
the basis of kinetic theory the power input from the bottom plate and the rate of the energy 
dissipation by inelastic collisions between particles. A balance between these two quantities 
at steady state determines a scaling of T on the system parameters. Since our theory can not 
estimate the energy dissipation by particle-particle collisions, it is incapable of deriving the 
dependence of c on r and N. 

We assumed that (r/r re i) 2 = (g/uoc) 2 <C 1, which assured that the system was in a 
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fluidized state, and Aquj/c <C 1, which allowed us to use hydrodynamic sound-wave theory. 
Since the acceleration F = Aqlu 2 jg > 1 in the fluidized state, we have g/uc < Aqlo/c, which 
implies that the second condition is dominant. 

We performed simulations with T > 10 and found that if Aqlo/c < 1, which corresponds 
to the case \[X < 1, the results agreed well with the theoretical predictions. They support 
our phenomenological argument deriving the Langevin equation (9), and our assumption that 
the granular fluid is well characterized by the total mass and the thermal velocity within the 
range of validity of our theory. 
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Fig. 3. (Color online) Comparison of the power spectrum Icm for the height of the center of mass 
obtained from simulations and the proposed theory for different N values: (a) N = 1000, (b) 
N = 100, (c) N = 10. The thick solid line depicts the theoretical predictions given by the second 
term in eq. (18) with fi — 2.0 and fl = 1.5. 
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